"""
    qgram - a program to browse and analyze earthquake data 
    using spectrograms 
    Copyright (C) 2011 John Ivey

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>

Contact the authors through the Google Code website http://code.google.com/p/qgram
"""
import numpy
import struct
from scipy.signal import cheby1, firwin, lfilter
import re

def read_psn(path):
    try:
        f = open(path,"rb")
    except:
        print("could not open ", path)
        raise
    # Identify type of file
    if f.read(10) == "PSNVOLUME1":
        # type = PSN volume file
        f.seek(0)
        bytes = f.read(12)
        volheader = struct.unpack("<10sh", bytes)
    else:
        f.seek(0)
    try:
        bytes = f.read(154)
        fixedheader = struct.unpack("<8slhcccccclddll3scccddccddd6s4s6sddhddd",bytes)
        filetype = fixedheader[0]
        variablelen = fixedheader[1]
        samplerate = fixedheader[11]
        datarecords = fixedheader[12]
        datatype = fixedheader[16]
        if filetype != "PSNTYPE4":
            print "file format", filetype, "not recognized"
            raise 
        f.read(variablelen)
        if datatype == '\x00':
            bytes1 = f.read(datarecords*2)
            formatstring = str(datarecords)+"h"
        elif datatype == '\x01':
            bytes1 = f.read(datarecords*4)
            formatstring = str(datarecords)+"i"
        elif datatype == '\x02':    
            bytes1 = f.read(datarecords*4)
            formatstring = str(datarecords)+"f"
        else:
            print("data type not recognized")
            raise
        try:
            data = struct.unpack(formatstring,bytes1)
        except:
            print("error unpacking data")
            raise
        floats = []
        for sample in data:
            floats.append(sample)
    except:
        print("error reading file")
        raise
                
    finally:
        f.close()
    
    x = numpy.array(floats)
    return x, samplerate

def read_text( path , skiprows = 0, column=0, delim = " " ):
    try:
        fread = open(path,"rb")
    except:
        print("could not open ", path)
        raise

    errors = 0
    nlines = 0
    ndata = 0
# the text file may be too big to read into memory so first count lines
    line = fread.readline()
    while line:
        if ( nlines >= skiprows ):
            ndata += 1
        line = fread.readline()
        nlines += 1

    cdata = 0
    linenum = 0
    
# preallocate numpy array - appending millions of times is slow
# array will be too big - will end up with extra zeros at end
# due to error lines - needs fixing
    output = numpy.zeros( ndata )
    fread.seek( 0 )

# skip past first lines in file    
    for i in range(0, skiprows):
        fread.readline()
        linenum += 1

# assume a float exists on each line, read and place in array         
    for clines in range( 0 , nlines ):
        line = fread.readline()
        # merge consecutive delimiters
        restring = delim + delim + "+"
        p = re.compile(restring)
        line1 = p.sub(delim, line)
        
        linepart = line1
#        for col in range(0,column+1):
#            valstr, dummy, linepart = linepart.partition(delim)
        
        
        try:
            valstr = line1.split(delim)[column]
            data = float(valstr)
            output[cdata] = data
            cdata += 1
        except:
            try:
                print "reached fallback delimited read code"
                data = float(linepart)
                cdata += 1
            except:
                errors += 1 
                print "error in line ",linenum
        linenum += 1 
        
    fread.close()
    if errors:
        print "errors ", errors
    print
    return output


def decimate_data(x, q, n=None, ftype='fir', axis=-1):
    """downsample the signal x by an integer factor q, using an order n filter
    
    By default, an order 8 Chebyshev type I filter is used or a 30 point FIR 
    filter with hamming window if ftype is 'fir'.

    (port to python of the GNU Octave function decimate.)

    Inputs:
        x -- the signal to be downsampled (N-dimensional array)
        q -- the downsampling factor
        n -- order of the filter (1 less than the length of the filter for a
             'fir' filter)
        ftype -- type of the filter; can be 'iir' or 'fir'
        axis -- the axis along which the filter should be applied
    
    Outputs:
        y -- the downsampled signal
    """

    if type(q) != type(1):
        raise "q should be an integer"

    if n is None:
        if ftype == 'fir':
            n = 30
        else:
            n = 8
    if ftype == 'fir':
        b = firwin(n+1, 1./(1.4*q), window='hanning')
        y = lfilter(b, 1., x, axis=axis)
    else:
        (b, a) = cheby1(n, 0.05, 0.8/q)
        y = lfilter(b, a, x, axis=axis)

    return y.swapaxes(0,axis)[n+1::q].swapaxes(0,axis)

